rm(list=ls())

#input data#
#setwd("/Volumes/MONOGAN/registrationSpatial/RRpolGeo/round2/")
laPanel<- read.csv("ccesLApanel.csv", header=T)

#load libraries
library(akima) #for drawing image/surface plots
library(geoR) #for variogram analyses
library(scatterplot3d) #for drop line plot

#preliminary OLS
year.fe<-lm(pid7~as.factor(year),data=laPanel)
summary(year.fe)

#define a "geodata" object
obj <- cbind(laPanel$eastings, laPanel$northings, laPanel$pid7)
#obj <- cbind(laPanel$eastings, laPanel$northings, year.fe$residuals)
la.geo <- as.geodata(obj, coords.col=1:2, data.col=3)

#Fit the exponential model using MLE
la.lik.fit <- likfit(la.geo, ini.cov.pars=c(6.0,1.0), cov.model = "exponential", trend = ~as.factor(laPanel$year), fix.nugget = FALSE, nugget = 1.0, nospatial = TRUE, lik.method = "ML") #note: "lik.method" is a correction of a typo from the text
summary(la.lik.fit)

#LR Test
chi2<-2*(-7239+7343);chi2#2(ln(full)-ln(null))
1-pchisq(chi2,df=2)

#compute correlations
exp(-1/6.578)
exp(-10/6.578)